Programmieren in R: Abgabe 1
Gabriel Torres Gamez, Si Ben Tran, Patrick Schürmann
Lineare und logistische Regressionen in R
Für die erste Abgabe im Modul PRR beschäftigen wir uns mit (multiplen) linearen und logistischen Regressionen. Diese Funktionen wenden wir auf den Datensatz “Carseats” an. Wir haben uns für diese Methoden entschieden, weil sie ein zentraler Bestandteil von Data Science sind. Sie wird verwendet, um Zusammenhängen in Datensätze zu verstehen und um Prognosen erstellen zu können. Diese Regressionen sind im R Grundpaket enthalten. Es muss also keine separate Library installiert und geladen werden.
Carseats ist ein generischer Datensatz, der im Buch “An Introduction to Statistical Learning” verwendet wird. Dieses Buch erarbeiten wir in diesem Studiengang im Modul “Statistisches Lernen”. Dort wird er zwar hauptsächlich angewendet, um Decision Trees zu erstellen. Dank 11 Variablen und 400 Observationen eignet er sich aber auch gut für Regressionen. Jede Zeile im Datensatz entspricht einem Verkaufslokal, gejoint mit Durchschnittsstatitistiken für das Einzugsgebiet des Lokals. Für eine genau Erklärung der Variablen verweisen wir auf die Homepage https://rdrr.io/cran/ISLR2/man/Carseats.html. Für die unten stehenden Hypothesen werden wir mit folgenden Variablen arbeiten:
- Sales: Verkaufte Anzahl Kindersitze pro Standort in 1’000 Stück
- Income: Das Lohnniveau der Einzugsgebiete der Standorte in 1’000 Dollar
- Advertising: Lokale Marketingausgaben in 1’000 Dollar
- Population: Bevölkerungsgrösse der Standorte in 1’000 Personen
- CompPrice: Die Preise für Kindersitze der lokalen Konkurrenten
- Price: Preis pro Kindersitz an den Standorten
- ShelveLoc: Ein Faktor der die Qualität des Produktes einteilt
- Age: Durchschnittsalter der lokalen Bevölkerung
- Education: Bildungslevel an jedem Standort
- Urban: Klassifikation, ob sich der Standort in einer städtischen oder ländlichen Region befindet
- US: Ein Indikator der angibt, ob sich das Geschäft in US befindet oder nicht
Lineare Regression
Bei der linearen Regression versucht man eine Zielvariable in Abhängigkeit durch eine andere Variable vorherzusagen. Das daraus entstehende Modell ist eine Gerade, welche durch folgende mathematische Formel beschrieben wird:
Formel: y = f(x) = ax + b = \(\beta_{0}\) + \(\beta_{1}\)x = mx + n
Das Wissen der linearen Regression, sowie der folgenden Regressionen, wird im Studiengang Data Science immersiv durch die Module llr, mlr und stl abgedeckt.
Lineare Multiple Regression
Bei der linearen multiplen Regression versucht man eine Zielvariable in Abhängigkeit durch mehrere andere Variablen vorherzusagen.
Formel ausgeschrieben: y = f(X) = \(\beta_{0}\) + \(\beta_{1}\) \(x_{1}\) + \(\beta_{2}\) \(x_{2}\) + … + \(\beta_{n}\) \(x_{n}\)
Formel vektorisiert: y = f(X) = X \(\cdot\) \(\vec{\beta_{n}}\) + \(\beta_{0}\)
Logistische Regression
Bei der logistischen Regression versucht man eine binäre Zielvariable in Abhängigkeit durch eine andere Variable vorherzusagen.
Sigmoid Formel: \(\sigma(z) = \frac{1}{1 + e^{-z}}\)
Formel fuer logistische Regression: f(x) = \(\sigma(ax + b)\) = \(\frac{1}{1 + e^{-(ax + b)}}\)
Auflistung der Hypothesen
Zur Analyse des Datensatzes haben wir uns folgende Fragen gestellt:
Sales ist abhängig von Income/Advertising/Population. (Lineare Regression)
Die Differenz von CompPrice und Price hat einen Einfluss auf den Verkauf. (Lineare Regression)
Sales lässt sich durch die Verwendung von mehreren Features zuverlässiger vorhersagen. (Multiple Lineare Regression)
Bewohner von urbanen Regionen kaufen weniger Kindersitze als Bürger ländlicher Regionen.
Vorbereitung
In diesem Kapitel installieren wir die notwendigen Libraries und verschaffen uns einen Überlick über den Carseats Datensatz.
Installationen
# Libraries installieren und laden
library(ISLR2)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.5 ✔ purrr 0.3.5
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()In der Library ISLR2 befindet sich der Datensatz. Tidyverse wird für die Darstellung des Datensatzes und Visualisierungen benötigt.
Überblick Datensatz Carseats
In einem ersten Schritt ist es wichtig den Datensatz zu verstehen, bevor Visualisierungen und Modelle erstellt werden.
Daten einlesen
# Carseats Daten einlesen
data(Carseats)Kopfzeile Carseats ausgeben
# Die ersten 6 Oberservationen ausgeben
head(Carseats)## Sales CompPrice Income Advertising Population Price ShelveLoc Age Education
## 1 9.50 138 73 11 276 120 Bad 42 17
## 2 11.22 111 48 16 260 83 Good 65 10
## 3 10.06 113 35 10 269 80 Medium 59 12
## 4 7.40 117 100 4 466 97 Medium 55 14
## 5 4.15 141 64 3 340 128 Bad 38 13
## 6 10.81 124 113 13 501 72 Bad 78 16
## Urban US
## 1 Yes Yes
## 2 Yes Yes
## 3 Yes Yes
## 4 Yes Yes
## 5 Yes No
## 6 No YesSelektieren von Spalten nach bestimmten Reihenfolgen
Carseats <- Carseats %>% select(Sales, Price, CompPrice, Advertising, ShelveLoc, Urban, US,
Population, Income, Age, Education,)Fusszeile des Datensatzes ausgeben
# Die letzten 6 Observationen ausgeben
tail(Carseats)## Sales Price CompPrice Advertising ShelveLoc Urban US Population Income Age
## 395 5.35 139 130 19 Bad Yes Yes 366 58 33
## 396 12.57 128 138 17 Good Yes Yes 203 108 33
## 397 6.14 120 139 3 Medium No Yes 37 23 55
## 398 7.41 159 162 12 Medium Yes Yes 368 26 40
## 399 5.94 95 100 7 Bad Yes Yes 284 79 50
## 400 9.71 120 134 0 Good Yes Yes 27 37 49
## Education
## 395 16
## 396 14
## 397 11
## 398 18
## 399 12
## 400 16Durch die Ausgabe von Kopfzeile und Fusszeile kommen wir mit dem Datensatz in Berührung und können uns ein Bild der Attribute sowie Observationen machen.
Informationen zum Datensatz herausfinden
# Informationen zum DataFrame ausgeben
str(Carseats)## 'data.frame': 400 obs. of 11 variables:
## $ Sales : num 9.5 11.22 10.06 7.4 4.15 ...
## $ Price : num 120 83 80 97 128 72 108 120 124 124 ...
## $ CompPrice : num 138 111 113 117 141 124 115 136 132 132 ...
## $ Advertising: num 11 16 10 4 3 13 0 15 0 0 ...
## $ ShelveLoc : Factor w/ 3 levels "Bad","Good","Medium": 1 2 3 3 1 1 3 2 3 3 ...
## $ Urban : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
## $ US : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...
## $ Population : num 276 260 269 466 340 501 45 425 108 131 ...
## $ Income : num 73 48 35 100 64 113 105 81 110 113 ...
## $ Age : num 42 65 59 55 38 78 71 67 76 76 ...
## $ Education : num 17 10 12 14 13 16 15 10 10 17 ...Durch str wird uns die Grösse des Dataframes angezeigt, sowie die Anzahl an Features (Spalten). Weiter erkennen wir von jeder Spalte den Datenyp und auch einige Werte, die in der Spalte vorkommen.
Reihenfolge von ShelveLoc definieren
Carseats$ShelveLoc <- factor(Carseats$ShelveLoc, ordered = TRUE, levels = c("Good", "Medium", "Bad"))Summarische Zusammenfassung von Carseats
# Statistische Kennzahlen
summary(Carseats)## Sales Price CompPrice Advertising ShelveLoc
## Min. : 0.000 Min. : 24.0 Min. : 77 Min. : 0.000 Good : 85
## 1st Qu.: 5.390 1st Qu.:100.0 1st Qu.:115 1st Qu.: 0.000 Medium:219
## Median : 7.490 Median :117.0 Median :125 Median : 5.000 Bad : 96
## Mean : 7.496 Mean :115.8 Mean :125 Mean : 6.635
## 3rd Qu.: 9.320 3rd Qu.:131.0 3rd Qu.:135 3rd Qu.:12.000
## Max. :16.270 Max. :191.0 Max. :175 Max. :29.000
## Urban US Population Income Age
## No :118 No :142 Min. : 10.0 Min. : 21.00 Min. :25.00
## Yes:282 Yes:258 1st Qu.:139.0 1st Qu.: 42.75 1st Qu.:39.75
## Median :272.0 Median : 69.00 Median :54.50
## Mean :264.8 Mean : 68.66 Mean :53.32
## 3rd Qu.:398.5 3rd Qu.: 91.00 3rd Qu.:66.00
## Max. :509.0 Max. :120.00 Max. :80.00
## Education
## Min. :10.0
## 1st Qu.:12.0
## Median :14.0
## Mean :13.9
## 3rd Qu.:16.0
## Max. :18.0Wir geben von jeder numerischen Spalte die summarischen Statistiken aus. Darin enthalten sind: Minimum, 1. Quantil, Median, Mittelwert, 3. Quantil, Maximum. Im summary() werden auch die Anzahl fehlende Werte angezeigt. In unserem Fall besitzt der Datensatz keine fehlende Werte. Die Aufbereitung der Daten entfällt deshalb.
Visualisierung der Daten
Verteilung der Sales
ggplot(data = Carseats, aes(x = Sales)) +
geom_density(fill = "lightgreen", color = "black", alpha = 0.3) +
labs(x = "Sales",
y = "Anzahl",
title = "Verteilung der Sales",
subtitle = "Carseats Datensatz")
In diesem Dichte-Diagramm erkennen wir die Verteilung der numerischen Variabel Sales. Die y Achse ist die Wahrscheinlichkeit für das Auftreten einer numerischen Sales Variabel. Aus diesem Grund wird die Dichtefunktion auch oft Wahrscheinlichkeitsdichtefunktion genannt. Sie ist eine geglättete Version des Histogramm und wird oft nach dem gleichen Konzept verwendet.
Verteilung der Sales anhand der Variabel Urban
ggplot(data = Carseats, aes(x = Sales, fill = Urban)) +
geom_density(alpha = 0.3) +
labs(x = "Sales",
y = "Anzahl",
title = "Verteilung der Sales",
subtitle = "Carseats Datensatz")
Analog zur ersten Visualisierung betrachten wir hier die Dichteverteilung von Sales mit der Auftrennung von Urban (ja/nein), sprich können wir mit diesem Plot Unterschiede zwischen urbanen Regionen (städtische oder ländliche) Sales betrachten. In unserem Fall gibt es keinen grossen Unterschied, man erkennt ganz leicht, dass Sales in nicht ländlichen Regionen diese im Sales Bereich grösser ist als bei städtischen Regionen.
Multiplots Verteilungen von Sales
ggplot(Carseats, aes(x = Sales, fill = US)) +
geom_density(alpha = 0.3) +
facet_grid(Urban ~ ShelveLoc, labeller = label_both) +
labs(x = "Sales",
y = "Anzahl",
title = "Verteilung der Sales unterteilt nach ShelveLoc und Urban",
subtitle = "Carseats Datensatz")
In diesem Multiplot unterteilen wir die Dichteverteilung nach weiteren Kategorien, wie Urban und ShelveLoc. Auch unterteilen wir farblich, ob sich das Geschäft in den US befindet oder nicht, damit wollen wir ein Gefühl für die Sales Verteilung mit verschiedenen kategorischen Werten kriegen. Aufällig ist, dass die Verteilungen bei ShelveLoc Bad deutlich weiter Links sind (Mittelwert tiefer als bei ShelveLoc Medium und Good). Auch erkennen wir bei ShelveLoc Good, dass die Kindersitze in den nicht US Ländern weniger verkauft werden als in den US selber.
Multiplots Verteilungen von Price
ggplot(Carseats, aes(x = Price, fill = US)) +
geom_density(alpha = 0.3) +
facet_grid(Urban ~ ShelveLoc, labeller = label_both) +
labs(x = "Price",
y = "Anzahl",
title = "Verteilung von Price unterteilt nach ShelveLoc und Urban",
subtitle = "Carseats Datensatz")
Analog zum Sales Multiplot erstellen wir ein Multiplot von Price, um eine Übersicht von Price mit unterschiedlichen Kategorien zu kriegen. Wir erkennen, dass bei ShelveLoc Good und urbanen Regionen die Preise in den US deutlich höher sind als in nicht US Ländern. Auch erkennen wir bei ShelveLoc Bad und nicht urbanen Regionen, das die Preise für einen Kindersitz bei nicht US Ländern teurer sind.
Hypothese 1: Sales ist abhängig von Income/Advertising/Population der Region.
Thema: Einfache lineare Regression
Sales vs. Income
Visualisierung von Sales und Income
plot(Sales ~ Income, data = Carseats,
main = "Sales vs Income", pch=21, bg = "lightblue")
Bevor wir mit der linearen Regression beginnen, erstellen wir zuerst ein Streudiagramm von beiden Variabeln, um die Zusammenhänge zu visualiseren. Hier erkennen wir im Plot, dass die Variabeln Einkommen und Sales nicht besonders stark korrelieren, da die Punkte überall im Plot verstreut sind.
Modell erstellen
lm.fit1 <- lm(Sales ~ Income, data = Carseats)Mit der Funktion “lm()” wird das Lineare Modell erstellt und trainiert.
Auswertung vom Modell
summary(lm.fit1)##
## Call:
## lm(formula = Sales ~ Income, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.2629 -1.9447 -0.1772 1.7654 8.9064
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.44356 0.37061 17.386 < 2e-16 ***
## Income 0.01533 0.00500 3.067 0.00231 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.795 on 398 degrees of freedom
## Multiple R-squared: 0.02309, Adjusted R-squared: 0.02063
## F-statistic: 9.407 on 1 and 398 DF, p-value: 0.00231R hat nun die lineare Regression erstellt. Die wichtigsten Werte lassen sich mit der Methode summary() anzeigen. Der Intercept beträgt 6.44 und die Steigung 0.015. Pro 1’000 Dollar höheres Einkommen steigt der Verkauf um 15 Stück. Dies erscheint im Verhältnis zum Intercept wenig. Es müsste darüber diskutiert werden, ob diese Erkenntnis für einen möglichen Kunden aussagekräftig ist. Anhand vom P Wert (0.00231) lässt sich aber feststellen, dass die Steigung siginfikant ist. Das Einkommen der Region hat also einen Einfluss auf die Verkaufszahlen. Vor einer endgültigen Aussage müssen aber weitere Untersuchungen, wie die Überprüfung der Residuen, durchgeführt werden. Diese Untersuchung führen wir später durch.
Übersicht der Modell Methoden
names(lm.fit1)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"Das Modell enthält weitere Berechnungen und Kennzahlen. Names() zeigt eine Übersicht dieser Methoden.
Koeffzienten vom Modell ausgeben
lm.fit1$coefficients## (Intercept) Income
## 6.44355748 0.01533361Die oben ersichtlichen Koeffizienten lassen sich direkt via $coefficients abrufen.
Konfidenzintervall für Koeffizienten angeben
confint(lm.fit1)## 2.5 % 97.5 %
## (Intercept) 5.714962797 7.17215216
## Income 0.005504876 0.02516235Wertvoll ist die Berechnung der Konfidenzintervalle, die im Summary nicht angezeigt wird. Diese gibt an, in welchem Bereich sich die Koeffizienten mit einer Wahrscheinlichkeit von 95% befinden.
Visualisierung der Daten mit dem Regressions Modell (Visualsierung ausbessern)
plot(Sales ~ Income, data = Carseats,
main = "Sales vs Income", pch=21, bg = "lightblue")
abline(lm.fit1, col = "red")
Das Modell lässt sich direkt mit der Plotfunktion darstellen. Mit abline() lassen sich Linien zeichnen, in diesem Fall die lineare Regression. Im Plot erkennen wir, dass Income und Sales eine Punktwolke bilden. Sales streut auf der ganzen X Achse von maximalen bis minimalen beobachteten Wert (0 - 15). Entsprechend korrelieren beide Variabeln wenig miteinander. Es gibt keinen klaren Zusammenhang.
Residuenanalyse
par(mfrow = c(2, 2))
plot(lm.fit1) 
Nun, da wir unser Modell haben, können wir feststellen, ob es ein gutes Modell ist, indem wir bei den Residuen folgende 3 Bedingungen überprüfen:
- Der Fehler sollte nicht abhängig sein
- Der Erwartungswert sollte bei 0 liegen.
- Die Residuen sollten einer Normalverteilung folgen
Das mittlere Kriterium lässt sich aus der Summary Funktion ablesen. Der Median der Residuen entspricht -0.18, was nahe genug ist. Dieses Kriterium ist erfüllt.
Auch dafür bietet R einfache Möglichkeiten an. Mit der folgenden beiden Funktionen lassen sich vier unterschiedliche Überprüfungen der Residuen anzeigen.
Im ersten Plot oben links (Tukeyanscombe Plot) wird ersichtich, dass die Residuen/Fehler unabhängig von einander sind und nicht keinem Muster folgen. Auch dieses Kriterium ist somit erfüllt.
Die Verteilung der Residuen lässt sich einfach im Histogramm überprüfen, ob deren Mittelwert sich um 0 verteilt.
Überprüfen vom Erwartungswert = 0 und Normalverteilt:
hist(lm.fit1$residuals,
main = "Verteilung der Residuen", col = "lightblue") 
Interpretation
Hypothese: Sales ist abhängig von Income
Wir erkennen aufgrund des Streudiagramm, dass eine Korrelation zwischen Einkommen und Sales kaum vorhanden ist. Das lineare Modell hat einen R^2 von 0.02, dies bestätigt die Aussage.
Aus diesem Grund ist Sales nicht abhängig von Income.
Sales vs. Population
Visualisierung von Sales und Population
plot(Sales ~ Population, data = Carseats,
main = "Sales vs Population", pch=21, bg = "lightblue")
Bevor das Modell erstellt wird, visualisieren wir beide Variabeln Population und Sales in einem Streudiagramm.
Model erstellen und Auswertung
lm.fit2 <- lm(Sales ~ Population, data = Carseats)
summary(lm.fit2)##
## Call:
## lm(formula = Sales ~ Population, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5864 -2.0176 -0.0597 1.6892 8.7213
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.2401837 0.2906658 24.909 <2e-16 ***
## Population 0.0009672 0.0009593 1.008 0.314
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.824 on 398 degrees of freedom
## Multiple R-squared: 0.002547, Adjusted R-squared: 4.116e-05
## F-statistic: 1.016 on 1 and 398 DF, p-value: 0.314Die Einfluss der Population ist nicht signifikant. Der folgende Plot bestätigt, dass die Steigung 0 beträgt.
Visualisierung der Daten mit dem Regressions Modell (Visualsierung ausbessern)
plot(Sales ~ Population, data = Carseats,
main = "Sales vs Population", pch=21, bg = "lightblue")
abline(lm.fit2, col = "red")
Hier in diesem Plot stellen wir die Population mit Sales sowie dem Modell zusammen. Wir erkennen, dass die Punkte sich deutlich von der roten Linie (Model) streuen. Das Problem besteht analog zu Sales und Income.
Residuenanalyse
par(mfrow = c(2, 2))
plot(lm.fit2) 
Auch in diesem Modell überprüfen wir die folgenden Kriterien des Models.
- Der Fehler sollte nicht abhängig sein
- Der Erwartungswert sollte bei 0 liegen.
- Die Residuen sollten einer Normalverteilung folgen
Das mittlere Kriterium lässt sich aus der Summary Funktion ablesen. Der Median der Residuen entspricht -0.0597, dieses Kriterium ist erfüllt.
Auch dafür bietet R einfache Möglichkeiten an. Mit der folgenden beiden Funktionen lassen sich vier unterschiedliche Überprüfungen der Residuen anzeigen.
Im ersten Plot oben links (Tukeyanscombe Plot) wird ersichtich, dass die Residuen/Fehler unabhängig von einander sind und nicht einem Muster folgen. Auch dieses Kriterium ist somit erfüllt.
Die Verteilung der Residuen lässt sich einfach im Histogramm überprüfen und ob sich deren Mittelwert um 0 verteilt.
hist(lm.fit2$residuals,
main = "Verteilung der Residuen", col = "lightblue") 
Interpretation
Hypothese: Sales ist abhängig von Population
Leider auch hier, ähnlich wie zu Sales vs. Income erkennen wir keinen Trend. Der R^2 Score beträgt 0.002547, es handelt sich also auch um kein verlässliches Modell.
Aus diesem Grund, stimmt auch diese Hypothese nicht.
Sales vs. Advertising
Visualisierung von Sales und Population
plot(Sales ~ Advertising, data = Carseats,
main = "Sales vs Advertising", pch=21, bg = "lightblue")
Bevor das Modell erstellt wird, visualisieren wir die beiden Variabeln Advertising und Sales.
Modell erstellen und Auswertung
lm.fit3 <- lm(Sales ~ Advertising, data = Carseats)
summary(lm.fit3)##
## Call:
## lm(formula = Sales ~ Advertising, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3770 -1.9634 -0.1037 1.7222 8.3208
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.7370 0.1925 35.007 < 2e-16 ***
## Advertising 0.1144 0.0205 5.583 4.38e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.723 on 398 degrees of freedom
## Multiple R-squared: 0.07263, Adjusted R-squared: 0.0703
## F-statistic: 31.17 on 1 and 398 DF, p-value: 4.378e-08Visualisierung der Daten mit dem Regressions Modell
plot(Sales ~ Advertising, data = Carseats,
main = "Sales vs Advertising", pch=21, bg = "lightblue")
abline(lm.fit3, col = "red")
Leider erkennen wir schon hier, dass es keinen Trend von Sales vs. Advertising vorhanden ist. Vollständigkeitshalber führen wir die Analysen trotzdem durch.
Residuenanalyse
par(mfrow = c(2, 2))
plot(lm.fit3) 
Auch in diesem Model überprüfen wir die folgenden Kriterien des Models.
- Der Fehler sollte nicht abhängig sein
- Der Erwartungswert sollte bei 0 liegen.
- Die Residuen sollten einer Normalverteilung folgen
Das mittlere Kriterium lässt sich aus der Summary Funktion ablesen. Der Median der Residuen entspricht -0.1037, dieses Kriterium ist erfüllt.
Auch dafür bietet R einfache Möglichkeiten an. Mit der folgenden beiden Funktionen lassen sich vier unterschiedliche Überprüfungen der Residuen anzeigen.
Im ersten Plot oben links (Tukeyanscombe Plot) wird ersichtich, dass die Residuen/Fehler unabhängig von einander sind und nicht keinem Muster folgen. Auch dieses Kriterium ist somit erfüllt.
Die Verteilung der Residuen lässt sich einfach im Histogramm überprüfen und ob sich deren Mittelwert um 0 verteilt.
hist(lm.fit3$residuals,
main = "Verteilung der Residuen", col = "lightblue")
Interpretation
Hypothese: Sales ist abhängig von Advertising
Leider ist auch bei Advertising und Sales kein linearer Zusammenhang erslichtlich. In der 3. Hypothese beschäftigen wir uns mit der multiplen linearen Regression, wir erwarten, dass sie besser sein wird.
Interpretation der Hypothese
Keines der linearen Modelle hat einen linearen Zusammenhang. Mit einem einfachen linearen Modell ist es nicht möglich, Sales vorherzusagen.
Hypothese 2: Die Differenz von CompPrice und Price hat einen Einfluss auf den Verkauf.
Thema: Einfach lineare Regression mit berechneten Daten
Sales vs. Preisdifferenz zu Konkurrenz
Berechnung der Spalte DiffPrice
Hinzufügen einer Spalte mit der Differenz zwischen unserem Preis und der Konkurrenz.
Carseats <- Carseats %>% mutate(DiffPrice = Price - CompPrice)
str(Carseats)## 'data.frame': 400 obs. of 12 variables:
## $ Sales : num 9.5 11.22 10.06 7.4 4.15 ...
## $ Price : num 120 83 80 97 128 72 108 120 124 124 ...
## $ CompPrice : num 138 111 113 117 141 124 115 136 132 132 ...
## $ Advertising: num 11 16 10 4 3 13 0 15 0 0 ...
## $ ShelveLoc : Ord.factor w/ 3 levels "Good"<"Medium"<..: 3 1 2 2 3 3 2 1 2 2 ...
## $ Urban : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 2 2 1 1 ...
## $ US : Factor w/ 2 levels "No","Yes": 2 2 2 2 1 2 1 2 1 2 ...
## $ Population : num 276 260 269 466 340 501 45 425 108 131 ...
## $ Income : num 73 48 35 100 64 113 105 81 110 113 ...
## $ Age : num 42 65 59 55 38 78 71 67 76 76 ...
## $ Education : num 17 10 12 14 13 16 15 10 10 17 ...
## $ DiffPrice : num -18 -28 -33 -20 -13 -52 -7 -16 -8 -8 ...Visualisierung von Sales und DiffPrice
plot(Sales ~ DiffPrice, data = Carseats,
main = "Sales vs DiffPrice", pch=21, bg = "lightblue")
In dieser Visualisierung von DiffPrice und Sales erkennen wir, dass die Streupunkte einem gewissen Trend folgen.
Model erstellen und Auswertung
lm.fit4 <- lm(Sales ~ DiffPrice, data = Carseats)
summary(lm.fit4)##
## Call:
## lm(formula = Sales ~ DiffPrice, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5885 -1.6583 -0.2678 1.5394 6.2367
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.691590 0.125567 53.29 <2e-16 ***
## DiffPrice -0.087662 0.005891 -14.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.267 on 398 degrees of freedom
## Multiple R-squared: 0.3575, Adjusted R-squared: 0.3559
## F-statistic: 221.5 on 1 and 398 DF, p-value: < 2.2e-16Die Auswertung zeigt, dass dieses Modell einen R^2 von 0.3575 hat. Dies ist nicht besonders gut, dennoch besser als die ersten drei Modelle.
Visualisierung der Daten mit dem Regressions Modell (Visualsierung ausbessern)
plot(Sales ~ DiffPrice, data = Carseats,
main = "Sales vs DiffPrice", pch=21, bg = "lightblue")
abline(lm.fit4, col = "red")
Hier in diesem Diagramm erkennen wir, dass das Modell eine negative Steigung hat. Dies erkennen wir einerseits an der Steigung im Summary, aber auch bei der Visualisierung.
Residuenanalyse
par(mfrow = c(2, 2))
plot(lm.fit4) 
Auch in diesem Modell uüberpruefen wir die folgenden Kriterien des Modells.
- Der Fehler sollten unabhängig sein.
- Der Erwartungswert sollte bei 0 liegen.
- Die Residuen sollten einer Normalverteilung folgen.
Das mittlere Kriterium lässt sich aus der Summary Funktion ablesen. Der Median der Residuen entspricht -0.2678, dieses Kriterium ist erfüllt.
Auch dafür bietet R einfache Möglichkeiten an. Mit der folgenden beiden Funktionen lassen sich vier unterschiedliche Überprüfungen der Residuen anzeigen.
Im ersten Plot oben links (Tukeyanscombe Plot) wird ersichtich, dass die Residuen/Fehler unabhängig von einander sind und nicht keinem Muster folgen. Auch dieses Kriterium ist somit erfüllt.
Die Verteilung der Residuen lässt sich einfach im Histogramm überprüfen und ob sich deren Mittelwert um 0 verteilt.
hist(lm.fit4$residuals,
main = "Verteilung der Residuen", col = "lightblue")
Interpretation der Hypothese
Hypothese: Die Differenz von CompPrice und Price hat einen Einfluss auf den Verkauf.
Die zweite Hypothese hat einen deutlich besseren R^2 Score verglichen zu den ersten drei Modellen bei der Hypothese 1. Die Differenz von CompPrice und Price hat aber keinen grossen Einfluss auf Sales, weil das Bestimmtheitsmass 0.3575 beträgt. Weil sich die Residuen also nur teilweise durch das Modell erklären lassen, verwerfen wir auch die zweite Hypothese.
Es ist interessant zu erkennen, dass der Preisunterschied zwischen der Konkurrenz und dem Produkt einen ganz leichten negativen Trend zeigt. Nichts destotrotz, sind wir mit diesem Modell nicht zufrieden und begeben uns direkt zur Hypothese 3.
Hypothese 3: Sales lässt sich durch die Verwendung von mehreren Features zuverlässiger vorhersagen.
Thema: Multiple lineare Regression
Sales vs. mehrere Features
Nachdem wir festgestellt haben, dass sich Sales nicht durch einzelne Variabeln erklären lasst, untersuchen wir die Möglichkeit, ob Sales sich durch mehrere Features erklären lässt. Wir werden sie zuerst gegenüber zwei Variabel untersuchen. Danach schauen wir uns den ganzen Datensatz an und behalten die aussagekräftigsten Variabeln.
Modelle erstellen und Auswertungen
In diesem Kapitel erstellen wir mehrere Modelle und wählen eines aufgrund einer ersten Auswertung für die weitere Untersuchung aus.
Sales vs. Income und Population
mlr.fit1 <- lm(Sales ~ Income + Population, data = Carseats)
summary(mlr.fit1)##
## Call:
## lm(formula = Sales ~ Income + Population, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.0371 -1.9342 -0.1477 1.7786 8.8532
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1785056 0.4493283 13.751 < 2e-16 ***
## Income 0.0153747 0.0049991 3.075 0.00225 **
## Population 0.0009902 0.0009493 1.043 0.29757
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.795 on 397 degrees of freedom
## Multiple R-squared: 0.02576, Adjusted R-squared: 0.02085
## F-statistic: 5.248 on 2 and 397 DF, p-value: 0.005627Zuerst untersuchen wir Sales gegenüber Income und der Population. In der Auswertung wird ersichtlich, dass Population keinen Zusammenhang aufweist. R^2 spricht ebenfalls gegen einen Zusammenhang in dieser Auswertung.
Die F-Statistik sagt aus, dass mit hoher Wahrscheinlichkeit (p-Wert < 0.6%) mindestens eine Variabel einen Einfluss hat. In unserem Fall ist das Income, was wir schon aus vorhergehenden Modellen wissen.
Sales vs. Income und Population und deren Multiplikation (Interaktion)
mlr.fit2 <- lm(Sales ~ Income * Population, data = Carseats)
summary(mlr.fit2)##
## Call:
## lm(formula = Sales ~ Income * Population, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9097 -1.9657 -0.2132 1.7682 8.8761
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.176e+00 7.950e-01 9.027 <2e-16 ***
## Income 1.161e-03 1.060e-02 0.110 0.913
## Population -2.655e-03 2.578e-03 -1.030 0.304
## Income:Population 5.198e-05 3.419e-05 1.520 0.129
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.79 on 396 degrees of freedom
## Multiple R-squared: 0.03141, Adjusted R-squared: 0.02408
## F-statistic: 4.281 on 3 and 396 DF, p-value: 0.005445Hier haben wir untersucht, ob die Variabeln Income und Population nicht nur additiv einen Einfluss haben, sondern ob sie auch gemeinsam wirken (Interaktion/Synergie). Das ist in diesem Fall aber nicht gegeben (P-Wert von 0.129 ist zu hoch).
Sales vs. alle Variabeln
mlr.fit3 <- lm(Sales ~ ., data = Carseats)
summary(mlr.fit3)##
## Call:
## lm(formula = Sales ~ ., data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8692 -0.6908 0.0211 0.6636 3.4115
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.9295889 0.5963763 13.296 < 2e-16 ***
## Price -0.0953579 0.0026711 -35.700 < 2e-16 ***
## CompPrice 0.0928153 0.0041477 22.378 < 2e-16 ***
## Advertising 0.1230951 0.0111237 11.066 < 2e-16 ***
## ShelveLoc.L -3.4295971 0.1082651 -31.678 < 2e-16 ***
## ShelveLoc.Q 0.3824279 0.0842896 4.537 7.62e-06 ***
## UrbanYes 0.1228864 0.1129761 1.088 0.277
## USYes -0.1840928 0.1498423 -1.229 0.220
## Population 0.0002079 0.0003705 0.561 0.575
## Income 0.0158028 0.0018451 8.565 2.58e-16 ***
## Age -0.0460452 0.0031817 -14.472 < 2e-16 ***
## Education -0.0211018 0.0197205 -1.070 0.285
## DiffPrice NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.019 on 388 degrees of freedom
## Multiple R-squared: 0.8734, Adjusted R-squared: 0.8698
## F-statistic: 243.4 on 11 and 388 DF, p-value: < 2.2e-16Mit dem Code “Sales ~ .” können wir ein Regressionsmodell Sales gegenüber allen Variabeln erstellen.
Dank der Auswertung sehen wir, dass sieben Variabeln in diesem Modell zum Sales beitragen. Population, Education und Location (Urban und US) haben keinen Einfluss. Wir werden nun die aussagekräftigsten vertieft überprüfen.
Die Fehlermeldung bei DiffPrice kommt daher, weil DiffPrice ein abhängies Feature von anderen Feature (Differenz aus CompPrice und Price) ist, welches wir in der vorherigen Hypothese erstellt haben. Das lineare Modell erkennt, dass DiffPrice in diesem Fall keine zusätzliche Information zum Modell beiträgt. Aus diesem Grund, weist uns die Auswertung ein NA aus.
Sales vs. aussagekräftigste Variabeln und DiffPrice
mlr.fit4 <- lm(Sales ~ Income + Advertising + ShelveLoc + Age + DiffPrice, data = Carseats)
summary(mlr.fit4)##
## Call:
## lm(formula = Sales ~ Income + Advertising + ShelveLoc + Age +
## DiffPrice, data = Carseats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7638 -0.6955 0.0266 0.6756 3.3990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.372761 0.222161 33.187 < 2e-16 ***
## Income 0.015903 0.001831 8.684 < 2e-16 ***
## Advertising 0.116005 0.007719 15.028 < 2e-16 ***
## ShelveLoc.L -3.415880 0.107703 -31.716 < 2e-16 ***
## ShelveLoc.Q 0.380870 0.083913 4.539 7.53e-06 ***
## Age -0.045853 0.003157 -14.522 < 2e-16 ***
## DiffPrice -0.095139 0.002660 -35.769 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.019 on 393 degrees of freedom
## Multiple R-squared: 0.8718, Adjusted R-squared: 0.8698
## F-statistic: 445.3 on 6 and 393 DF, p-value: < 2.2e-16Durch die Auswahl der aussagekräftigsten Variabeln konnten wir ein Modell erstellen, das Sales zuverlässig erklärt. Die verwendeten Variabeln haben wir aus dem Modell entnommen, welche einen tiefen P-Wert aufweisen. In diesem Modell haben wir nur (auch in der F-Statistik) sehr tiefe p-Werte und ein ziemlich hohes R^2 von 0.87. Wir überprüfen deshalb nun die Residuen.
Residuenanalyse
par(mfrow = c(2,2))
plot(mlr.fit4)
Im Q-Q Plot wird bestätigt, dass es sich um ein lineares Modell handelt. Der erste Plot und das folgende Histogram bestätigen ebenfalls, dass die Residuen normalverteilt und einen Mittelwert von 0 aufweisen.
hist(mlr.fit4$residuals,
main = "Verteilung der Residuen", col = "lightblue")
Interpretaion der Hypothese
Hypothese: Sales lässt sich durch die Verwendung von mehreren Features zuverlässiger vorhersagen.
Unser letztes Modell hat gezeigt, dass mit einer sorgfältigen Auswahl der Variabeln, Sales gut und zuverlässig erklärt werden kann. In der Praxis könnten noch mehr Kombinationen (Variabeln, Interaktionen, exponentielle Zusammenhänge etc…) untersucht werden. In dieser Abgabe lassen wir das weg, weil es sich um Wiederholungen der gezeigten Funktionen handelt.
Hypothese 4: Bewohner von urbanen Regionen kaufen weniger Kindersitze als Bürger ländlicher Regionen.
Thema: Logistische Regression
Sales vs. Urban
Klassifizierung erstellen
Um logistische Regressionen durchzuführen, benötigen wir binäre Variabeln. Im ersten Schritt ersetzen wir Yes/No von Urban durch 1/0.
Carseats <- Carseats %>%
mutate(Urban_binaer = recode(Urban, "No" = 0, "Yes" = 1))Visualisierung von Sales und Urban
ggplot(Carseats, aes(x=Sales, y=Urban)) +
geom_point(color = "blue") +
labs(x = "Sales",
y = "Urban",
title = "Sales vs Urban",
subtitle = "Carseats Datensatz")
Da wir in den Plots zur logistischen Regression die Sigmoid Funktion einzeichnen möchten, verwenden wir diese Darstellung bei sämtliche logistischen Regression.
Dank der visuellen Analyse von Sales und Urban wird ersichtlich, dass Urban:Yes in beide Richtungen stärker streut als Urban:No. Entsprechend wird das Modell kaum verlässliche Werte ergeben.
Model erstellen und Auswertung
logreg.fit1 <- glm(Urban_binaer ~ Sales, data = Carseats , family = binomial)
summary(logreg.fit1)##
## Call:
## glm(formula = Urban_binaer ~ Sales, family = binomial, data = Carseats)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6032 -1.5464 0.8288 0.8396 0.8737
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.96122 0.31250 3.076 0.0021 **
## Sales -0.01197 0.03883 -0.308 0.7578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 485.25 on 399 degrees of freedom
## Residual deviance: 485.16 on 398 degrees of freedom
## AIC: 489.16
##
## Number of Fisher Scoring iterations: 4Wie wir aufgrund der Visualisierung erkannt haben, ergibt ein logistisches Regressionsmodell einen grossen P>|z| Wert. Somit gibt es keinen logistischen Zusammenhang zwischen Sales und Urban.
Gegenüber der lineare Regression haben wir den Parameter “family = binomial” erfasst.
Sales vs. ShelveLoc
Um ein Beispiel für eine logistische Regression zeigen zu können, untersuchen wir auch Sales und ShelveLoc
Visualisierung von Sales und ShelveLoc
Zuerst für alles ShelveLoc, danach haben die Location “Medium” entfernt.
Carseats %>%
ggplot(aes(x = Sales, y = ShelveLoc)) +
geom_point(color = "blue") +
labs(x = "Sales",
y = "ShelveLoc",
title = "Sales vs ShelveLoc",
subtitle = "Carseats Datensatz")
Carseats_ex_medium <- Carseats %>% filter(ShelveLoc != "Medium")
Carseats_ex_medium %>%
ggplot(aes(x = Sales, y = ShelveLoc)) +
geom_point(color = "blue") +
labs(x = "Sales",
y = "ShelveLoc",
title = "Sales vs ShelveLoc",
subtitle = "Carseats Datensatz")
Im ersten Plot wurden drei Variabeln dargestellt. Für ein gutes Vorzeigebeispiel verwenden wir nur Good und Bad.
Klassifizierung erstellen
Carseats_ex_medium <- Carseats_ex_medium %>%
mutate(ShelveLoc_binaer = recode(ShelveLoc, "Bad" = 0, "Good" = 1))
head(Carseats_ex_medium)## Sales Price CompPrice Advertising ShelveLoc Urban US Population Income Age
## 1 9.50 120 138 11 Bad Yes Yes 276 73 42
## 2 11.22 83 111 16 Good Yes Yes 260 48 65
## 5 4.15 128 141 3 Bad Yes No 340 64 38
## 6 10.81 72 124 13 Bad No Yes 501 113 78
## 8 11.85 120 136 15 Good Yes Yes 425 81 67
## 11 9.01 100 121 9 Bad No Yes 150 78 26
## Education DiffPrice Urban_binaer ShelveLoc_binaer
## 1 17 -18 1 0
## 2 10 -28 1 1
## 5 13 -13 1 0
## 6 16 -52 0 0
## 8 10 -16 1 1
## 11 10 -21 0 0Wir haben “Medium” weggelassen und wieder eine binäre Klassifizierung erstellt.
Model erstellen und Auswertung
logreg.fit2 <- glm(ShelveLoc_binaer ~ Sales, data = Carseats_ex_medium , family = binomial)
summary(logreg.fit2)##
## Call:
## glm(formula = ShelveLoc_binaer ~ Sales, family = binomial, data = Carseats_ex_medium)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3930 -0.5021 -0.1506 0.5015 2.5892
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.0248 0.8714 -6.914 4.71e-12 ***
## Sales 0.7566 0.1082 6.990 2.76e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 250.25 on 180 degrees of freedom
## Residual deviance: 136.49 on 179 degrees of freedom
## AIC: 140.49
##
## Number of Fisher Scoring iterations: 5Ein tiefer P-Wert weist auf eine hohe Zuverlässigkeit vom Modell hin.
Im folgenden Plot wird ersichtlich, dass das Modell bei allen Läden mit Sales > 7.96 den ShelveLoc als Good klassifiziert. Die Berechnung des Wertes haben wir im folgenden Code Chunk auskommentiert.
Visualisierung der Daten mit dem Regressions Modell
ggplot(Carseats_ex_medium, aes(x=Sales, y=ShelveLoc_binaer)) +
geom_point(color = "blue") +
labs(x = "Sales",
y = "ShelveLoc",
title = "Sales vs ShelveLoc",
subtitle = "Carseats Datensatz") +
stat_smooth(method="glm", formula = y ~ x, color="red",
se = FALSE, method.args = list(family=binomial)) +
geom_vline(xintercept = 7.9633)
# Berechnung der Decision Boundary
#-logreg.fit2$coefficients[1] / logreg.fit2$coefficients[2]Evaluierung vom logistischen Model mittels Confusions Matrix
Um die Metriken des Modelles zu berechnen, erstellen wir eine Confusion Matrix.
Confusion Matrix erstellen
# zuerst erstellen wir eine Prediction, mit den Testdaten
# Predicte ShelveLoc_binaer Werte
logreg.fit2.probs <- predict(logreg.fit2, type = "response")
# Erstellt einen Vektor mit 181 Null Werten
logreg.fit2.pred <- rep(0, 181)
# Uebreschreibt alle Werte die groesser sind als 0.5 als 1
logreg.fit2.pred[logreg.fit2.probs > .5] = 1
# Erstellen eine Kreuztabelle (Confusions Matrix)
table("Actual Class" = Carseats_ex_medium$ShelveLoc_binaer, "Predicted Class" = logreg.fit2.pred)## Predicted Class
## Actual Class 0 1
## 0 82 14
## 1 19 66Die Spalten entsprechen den tatsächlichen Werten und in den Zeilen befinden sich die vorhergesagten Werte. Lesebeispiel: Das Modell hat 19 Werte als 0/Bad klassifiziert, die aber tatsächlich ein 1/Good gewesen sind.
Confusion Matrix Performenz Kennzahlen
Accuracy
Die Accuracy entspricht dem Anteil korrekter Prognosen an allen Prognosen.
# Accuracy Berechnen
acc <- (82 + 66) / (82 + 19 + 14 + 66)
acc## [1] 0.8176796Precision
Precision = True Positiv / (True Positiv + False Positiv) Preicision ist, Anteil richtig positiv an allen Positiv
Precision entspricht der Anteil korrekter positiven/1 Prognosen an allen positiven Prognosen.
# Precision Berechnen
precision <- 66 / (66 + 14)
precision## [1] 0.825Recall / Sensitivity
Recall = True Positiv / (True Positiv + False Negativ) Recall ist, Anteil richtig positiv an allen richtigen Vorhersagen
Recall / Sensitivity entspricht der Anteil korrekter positiven/1 Prognosen an allen tatsächlich positiven Werten.
# Recall berechnen
recall <- 66 / (66 + 19)
recall## [1] 0.7764706Zusammenfassung
Mit unseren verwendeten Modellen konnten wir aufzeigen, dass Sales nicht von einzelnen Variabeln, aber von mehreren zusammenhängt. Dank der multiplen linearen Regression konnten wir zuverlässige Einflussfaktoren aufzeigen.
Dank der logistischen Regression konnten wir darstellen, dass Urban keinen Einfluss hat. Mit dem Weglassen von Medium konnten wir für ShelveLoc ein zuverlässiges Modell erstellen.